clear; close all; clc; addpath('../utilities');

para = [0.043682298946	10.680407044771	0.031532285160	25.001413406635	3.812243088585	0.203429573303	0.102423908057	7.945920271541];

main_para

options = optimoptions('fsolve','Display','off','TolFun',1e-12);
beta    = fsolve(@solve_beta,0.99,options,gamma,psi,mu_C,sig_C,P,P_nodis);

main_SDF

lamS    = ones(nS,nK);
lamD    = ones(nS,nK);
k       = linspace(k_min,k_max,nK)';
k_bar   = ones(nS,1)*5;
k_iss   = ones(nS,1)*10;

[kgrid,mgrid] = meshgrid(k,(1:nS)');
kgrid         = kgrid(:);
mgrid         = mgrid(:);

Q_big   = repmat(Q,nK,1);
Rf1_mat = repmat(Rf(:,1),1,nK);

disp('Solve VFI')
tic, main_VFI, minutes = toc/60

disp('Solve CDS pricing')
main_CDS
main_CDS_P

disp('Solve option pricing')
money       = ones(nS,nK);
main_options

IV          = NaN(nS,nK,2);
IV(:,:,1)   = BSIV;

money       = exp(-BSIV/sqrt(12));
main_options

IV(:,:,2)   = BSIV;
IV_ATM      = IV(:,:,1)';
IV_ATM      = permute(IV_ATM,[2,1]);
IV_Skew     = IV(:,:,2)'-IV(:,:,1)';
IV_Skew     = permute(IV_Skew,[2,1]);

disp('Solve cond. return moments')
main_Ret

disp('Solve cond. return moments')
main_options_on_assets
save model_solution omega stat stat_Q ATM_Asset_IV

disp('Run simulation')
tic, main_simu, minutes = toc/60
save model_simulation table_5 table_7 table_8 ts_var cs_var states G_m

disp('Run Figure 7 calculations')
tic, Figure_7_calculations, minutes = toc/60
save Figure_7_out EIV Asset_IV Asset_IV_noIdio KS_vec MON_values prob_sim

disp('Run ergodic simulation')
tic, main_simu_pop, minutes = toc/60
save model_simulation_pop table_5_pop table_7_pop

disp('Run crisis simulation for Tables 3-5')
tic, main_simu_crisis, minutes = toc/60
save crisis_simulation_out MOM1 MOM2 MOM3 MOM4 MOM4 MOM5 MOM6 MOM7 MOM8 MOM9 MOM10 MOM11 state1 state2 Dsize Dmin Dmax para  

disp('Run Figure 6 simulation')
tic, Figure_6_simulation, minutes = toc/60
save Figure_6_out LEV VOL CDS cut




